source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

inflam15 <- read_csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

modc_marker <- list('cMono'=c('CD14','S100A8','S100A9'),
                    'ncMono' = c('FCGR3A','CDKN1C'),
                    'Macrophage'=c('CD68','MAFB','MARCO','CD163','MRC1'),
                    'cDC'=c('FCER1A','CLEC9A','XCR1','CD1C','CLEC10A'),
                    'pDC'='TLR7')

# 3v3 data ----------
sobj <- read_rds('mission/FPP/psoriasis/gao2021harbin.rds')

sobj <- sobj |>
  mutate(hpca_main = case_when(hpca_main == 'CMP' ~ 'Mast_cells',
                               hpca_main == 'Monocyte' ~ 'Macrophages',
                               .default = hpca_main))

sobj |> FindMarkers(ident.1 = 'CMP') |>
  as_tibble(rownames = 'gene')

sobj |>
  DimPlot(group.by = 'hpca_main', cols = 'Paired') +
  ggtitle('Psoriasis skin')

sobj |>
  DotPlot(modc_marker, cols = 'RdBu')

sobj |>
  DotPlot2d('TRPM2', group, hpca_main) +
  labs(x = 'Group', y = 'Cell type')

sobj_m <- sobj |> filter(hpca_main == 'Monocyte') |>
  quick_process_seurat(skip_norm = T)

sobj_m |> FeaturePlot('TRPM2', split.by = 'group', order = T,
                      cols = c('lightgrey','red'))

sobj_m |>
  DotPlot(list_c(modc_marker), cols = 'RdBu') +
  RotatedAxis()

sobj_m <- sobj_m |>
  mutate(manual_fine = ifelse(seurat_clusters %in% c(2,5,6,7),
                              'DC', 'Macrophages'))

sobj_m |> DimPlot(group.by = 'manual_fine') +
  ggtitle('DC/Macrophages')

psovhc_myl_deg <- sobj_m |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'group',
                       ident.1 = 'Psor')

psovhc_myl_deg |>
  filter(gene == 'TRPM2')

m2_bysample <- sobj_m |>
  DotPlot2d('TRPM2', manual_fine, orig.ident) |>
  pluck('data') |>
  as_tibble()
  
m2_bysample |>
  mutate(group = str_remove(group.y, '\\d')) |>
  filter(group.x != 'DC') |>
  ggplot(aes(group, avg.exp, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(method = 't.test', size = 2,
                     comparisons = list(c('Ctrl','Psor'))) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  scale_fill_hue(direction = -1) +
  labs(title = 'TRPM2 in skin macrophage', y = 'Average expression')

publish_pdf('gao21.psor.mf.pdf')

sobj_m |>
  filter(manual_fine != 'DC') |>
  bill.violin('TRPM2', group.by = group) +
  labs(y = 'Normalized expression', title = 'TRPM2 in skin macrophages') +
  scale_fill_hue(direction = -1)

sobj_mf <- sobj_m |>
  filter(manual_fine != 'DC') |>
  FindClusters(algorithm = 4, random.seed = 1)

sobj_mf |>
  DotPlot('TRPM2', cols = 'RdBu')

sobj_mf <- sobj_mf |>
  mutate(trpm2_type = ifelse(seurat_clusters %in% c(1,4,5),
                             'TRPM2-hi MF', 'TRPM2-lo MF'))

sobj_mf |>
  filter(trpm2_type != 'TRPM2-hi MF') |>
  bill.violin('TRPM2', group.by = group) +
  labs(y = 'Normalized expression', title = 'TRPM2 in TRPM2-high macrophages') +
  scale_fill_hue(direction = -1)

# GSE173706 ----------
shelf(GEOquery)

getGEOSuppFiles('GSE173706', makeDirectory = F, fetch_files = F)

merl_path <-
  list.files('mission/SLE_TRPM2_MfMo/data/merleev22psoriasis/', full.names = T)

import_ensembl_mat <- function(path) {
  foo <- data.table::fread(path)
  
  orgeg <- foo$V1 |>
    clusterProfiler::bitr(fromType = 'ENSEMBL', toType = 'SYMBOL',
                          OrgDb = 'org.Hs.eg.db')

  foo |>
    right_join(orgeg, join_by(V1 == ENSEMBL)) |>
    distinct(SYMBOL, .keep_all = T) |>
    mutate(V1 = NULL) |>
    column_to_rownames('SYMBOL') |>
    as.sparse()
}

merl_lst <- merl_path |>
  map(import_ensembl_mat, .progress = T)

merl_lst |> str()

merl_name <- merl_path |>
  str_extract('GSM.+(?=.csv)') |>
  str_remove('GSM\\d+_')

merl_lst <- merl_lst |>
  map2(merl_name, add_name_suffix)

new.mat <- merl_lst[[1]] |>
  RowMergeSparseMatrices(merl_lst[-1])

sobj <- new.mat |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

merl_name

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'NS|PP|PN')) |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito_ratio < 10) |>
  quick_process_seurat()

sobj |>
  write_rds('mission/SLE_TRPM2_MfMo/merl22psoria.rds')

sobj <-
  read_rds('mission/SLE_TRPM2_MfMo/merl22psoria.rds')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca_main')

sobj |> DotPlot(pbmc_markers, cols = 'RdBu')

sobj <- sobj |>
  mutate(manual_main = case_when(hpca_main == 'DC' ~ 'Macro/DC',
                                 hpca_main == 'CMP' ~ 'Mast_cells',
                                 seurat_clusters == 13 ~ 'Melanocytes',
                                 str_detect(hpca_main, 'stem') ~ 'Fibroblasts',
                                 .default = hpca_main))

## total umap ------------
sobj |>
  ggplot(aes(umap_1, umap_2, fill = manual_fine)) +
  geom_bin2d(bins = 512) +
  labs(title = 'Psoriasis skin', fill = '') +
  scale_fill_brewer(palette = 'Paired') +
  theme_jpub(theme_classic)

publish_pdf('psoria.skin.umap.pdf', width = 70)

sobj_myl <- sobj |>
  filter(manual_main == 'Macro/DC')

sobj_myl <- sobj_myl |>
  quick_process_seurat(skip_norm = T)

sobj_myl |>
  DotPlot(c(modc_marker, 'CD207'), cols = 'RdBu', cluster.idents = T) +
  RotatedAxis()

sobj_myl <- sobj_myl |>
  filter(!(seurat_clusters %in% c(12,8,10)))

sobj_myl <- sobj_myl |>
  FindClusters(algorithm = 4, random.seed = 1)

sobj_myl |> DimPlot(cols = 'Paired', label = T, label.box = T,
                    repel = T, label.size = 2) +
  theme_jpub(theme_classic)

publish_pdf('psoria.skin.myeloid.umap.pdf', width = 70)

sobj_myl <- sobj_myl |>
  mutate(manual_fine = case_when(seurat_clusters %in% c(5,9,2,4) ~ 'Macro',
                                 seurat_clusters %in% c(1,6) ~ 'LC',
                                 .default = 'DC'))

sobj_myl |> DimPlot(group.by = 'manual_fine', cols = 'Paired') +
  ggtitle('Skin DC/macrophages')

sobj_myl |>
  DotPlot(c('TRPM2')) |>
  pluck('data') |>
  BubblePlot() +
  theme_jpub()

publish_source_plot('psoria.skin.myeloid.m2.dotplot')

## macro ---------
sobj_mf <- sobj_myl |>
  filter(manual_fine == 'Macro') |>
  quick_process_seurat(skip_norm = T)

sobj_mf <- sobj_mf |>
  AddModuleScore(features = list(inflam15$gene), name = 'inflam')

sobj_mf |> write_rds('mission/SLE_TRPM2_MfMo/merl22psoria_macro.rds')

sobj_mf <- read_rds('mission/SLE_TRPM2_MfMo/merl22psoria_macro.rds')

sobj_mf |>
  DotPlot(c('TRPM2')) |>
  pluck('data') |>
  BubblePlot() +
  theme_jpub()

publish_source_plot('psoria.skin.mf.m2.dotplot')

inflam_m2 <- sobj_mf |>
  DotPlot(c('inflam1', 'TRPM2', 'CCR2', 'HMGB1')) |>
  pluck('data') |>
  as_tibble()

inflam_m2 |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .12, 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2, label.y = .9) +
  labs(title = 'Macrophage clusters in psoriasis skin',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('psoria.m2.inflam.corr', width = 70)

## psoria vs hc ---------
hcvo_mf_deg <- sobj_mf |>
  FindMarkers(group.by = 'group', ident.1 = 'NS')

hcvo_mf_deg |>
  as_tibble(rownames = 'gene') |>
  filter(gene == 'TRPM2')

sobj_mf |>
  bill.violin('TRPM2', group)

sobj_myl |>
  DotPlot2d('TRPM2', group, manual_fine)

mf_m2_bysample <- sobj_mf |>
  DotPlot('TRPM2', group.by = 'orig.ident') |>
  pluck('data')

mf_m2_bysample |>
  mutate(group = str_extract(id, 'PP|PN|NS') |>
           case_match('NS' ~ 'HC', 'PP' ~ 'Psoriasis', .default = 'PN')) |>
  filter(group != 'PN') |>
  ggplot(aes(group, avg.exp, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  scale_fill_hue(direction = -1) +
  labs(title = 'TRPM2 in skin macrophage', y = 'Average expression') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_compare_means(comparisons = list(c('HC','Psoriasis')),
                     method = 't.test', size = 2)

publish_source_plot('psoriasis.mer22.mf.trpm2.bysample')

sobj_mf |>
  DimPlot()

sobj_mf |>
  filter(group != 'PN') |>
  mutate(group = case_match(group,
                            'NS' ~ 'HC', 'PP' ~ 'Psoriasis', .default = 'PN')) |>
  FeaturePlot('TRPM2', split.by = 'group', order = T,
              cols = c('lightgrey','red')) &
  theme_jpub(theme_classic) & NoLegend()

publish_pdf('psoria.skin.mf.m2.featureplot.pdf', width = 100)

mf_frac_in_skin <- sobj |>
  mutate(manual_main = ifelse(.cell %in% colnames(sobj_mf), 'MF', 'other')) |>
  calc_frac_conf_on_grouped_count(orig.ident, manual_main)

mf_frac_in_skin |>
  mutate(group = str_extract(orig.ident, 'PP|PN|NS') |>
           case_match('NS' ~ 'HC', 'PP' ~ 'Psoriasis', .default = 'PN')) |>
  filter(manual_main == 'MF', group != 'PN') |>
  ggplot(aes(group, fraction*100, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  scale_fill_hue(direction = -1) +
  labs(title = 'Macrophage fraction in skin', y = '% in total cells') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_compare_means(comparisons = list(c('HC','Psoriasis')),
                     method = 't.test', size = 2)

publish_source_plot('psoriasis.mer22.mf.frac.bysample')

sobj <- sobj |>
  mutate(manual_main = case_when(.cell %in% colnames(sobj_mf) ~ 'Macrophage',
                                 manual_main == 'Macro/DC' ~ 'DC',
                                 .default = manual_main))

sobj |>
  mutate(group = str_extract(orig.ident, 'PP|PN|NS') |>
           case_match('NS' ~ 'HC', 'PP' ~ 'Psoriasis', .default = 'PN')) |>
  filter(group != 'PN') |>
  DotPlot2d('TRPM2', group, manual_main)

psor_m2_dotplot <- last_plot() |>
  pluck('data')

psor_m2_dotplot |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 in skin') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('psoria.merl22.skin.m2.dotplot', width = 55)

sobj_mf |>
  DotPlot(c('TRPM2','CD14'), cols = 'RdBu')

m2h_bc <- sobj_mf |>
  filter(seurat_clusters == 4) |>
  colnames()

m2l_bc <- sobj_mf |>
  filter(seurat_clusters != 4) |>
  colnames()

sobj <- sobj |>
  mutate(manual_fine = case_when(.cell %in% m2h_bc ~ 'CD14+ Macrophage',
                                 .cell %in% m2l_bc ~ 'CD14- Macrophage',
                                 manual_main == 'Macro/DC' ~ 'DC',
                                 .default = manual_main))

sobj |>
  mutate(group = str_extract(orig.ident, 'PP|PN|NS') |>
           case_match('NS' ~ 'HC', 'PP' ~ 'Psoriasis', .default = 'PN')) |>
  filter(group != 'PN') |>
  DotPlot2d('TRPM2', group, manual_fine)

psor_m2_dotplot <- last_plot() |>
  pluck('data')

psor_m2_dotplot |>
  mutate(group.y = fct_reorder(group.y, avg.exp, max)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 in skin') +
  theme_jpub()

publish_source_plot('psoria.cd14mf.skin.m2.dotplot', width = 60)

## M2-hi vs M2-lo --------
sobj_mf |>
  filter(group == 'PP') |>
  mutate(trpm2_type = ifelse(seurat_clusters == 1, 'TRPM2-hi', 'TRPM2-lo')) |>
  bill.violin(inflam15$gene, trpm2_type, facet.ncol = 5) +
  theme(axis.text.x = element_blank())

# PBMC GSE194315 -------
psa_meta <- 
  read_delim('mission/SLE_TRPM2_MfMo/data/liu22psa/GSE194315_CellMetadata-PSA_TotalCiteseq_20220103.tsv.gz')

psa_meta$Status |> unique()

supp_name <- list.files('mission/SLE_TRPM2_MfMo/data/liu22psa/', 'PBMC',
                        full.names = T)

# ~7min for 2 worker read in
psa_mex <- read_geo_supp(supp_name, name_regex = 'PBMC.{5}', workers = 2)

psa_mex

meta_bc <- psa_meta |>
  filter(IncludedInStudy) |>
  pull(CellName) |>
  map_chr(\(x)str_c(x, '-1'))

psa_mex <- psa_mex |>
  filter(.cell %in% meta_bc)

sobj <- psa_meta |>
  filter(IncludedInStudy) |>
  mutate(.cell = str_c(CellName, '-1'), .keep = 'unused') |>
  left_join(x = psa_mex, y = _)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito_ratio < 10)

# already annotated, saved umap embedding
sobj <- sobj |>
  NormalizeData()

sobj <- sobj |>
  quick_process_seurat()

sobj <- sobj |>
  mutate(manual_main = case_when(str_detect(CellType, 'Mono') ~ 'Monocytes',
                                 str_detect(CellType, 'cDC') ~ 'cDC',
                                 str_detect(CellType, 'DC') ~ 'pDC',
                                 str_detect(CellType, 'NK') ~ 'NK',
                                 str_detect(CellType, '^B|Plasma') ~ 'B cells',
                                 .default = 'T cells',))

sobj <- sobj |>
  filter(!(CellType %in% c('Doublet','Eryth','Platelet')))

sobj |>
  ggplot(aes(x = umap_1, y = umap_2, fill = manual_main)) +
  geom_bin2d(bins = 512) +
  scale_fill_brewer(palette = 'Paired') +
  labs(fill = 'Cell type', title = 'Psoriasis patient PBMC') +
  theme_jpub(theme_classic)

sobj |> ggplot(aes(x = umap_1, y = umap_2, fill = CellType)) +
  geom_bin2d(bins = 64) + theme_classic() +
  scale_fill_manual(values = DiscretePalette(36))

sobj <- sobj |>
  filter(Status != 'PSX') |>
  mutate(manual_coarse = ifelse(str_detect(manual_main, 'DC'), 'DC', manual_main),
         Status = case_match(Status, 'PSA' ~ 'Psoriasis+Arthritis',
                             'PSO' ~ 'Only psoriasis', .default = Status),
         Disease = ifelse(Status == 'Healthy', 'Healthy', 'Psoriasis'))

sobj |>
  as_tibble() |>
  write_source_csv('psoria_pbmc_meta')

pso_pbmc_meta <-
  read_csv('mission/SLE_TRPM2_MfMo/results/psoria_pbmc_meta.csv')

pso_pbmc_meta |>
  filter(CellType != 'pDC') |>
  ggplot(aes(x = umap_1, y = umap_2, fill = manual_coarse)) +
  geom_bin2d(bins = 512) +
  theme_jpub(theme_classic) +
  scale_fill_brewer(palette = 'Paired') +
  labs(fill = 'Cell type', title = 'PBMC')

pso_pbmc_meta |>
  filter(CellType != 'pDC') |>
  mutate(manual_fine = case_when(.cell %in% m2l_bc ~ 'CD16+ Monocyte',
                                 manual_main == 'Monocytes' ~ 'CD14+ Monocyte',
                                 .default = manual_coarse)) |>
  ggplot(aes(x = umap_1, y = umap_2, fill = manual_fine)) +
  geom_bin2d(bins = 512) +
  theme_jpub(theme_classic) +
  scale_fill_brewer(palette = 'Paired') +
  labs(fill = 'Cell type', title = 'PBMC')

publish_pdf('psoria.pbmc.umap.pdf', width = 70)

## TRPM2 in pbmc ----------
sobj |>
  DotPlot2d('TRPM2', Status, manual_main) +
  labs(x = 'Group', y = 'Cell type') +
  theme_jpub() +
  scale_radius(range = c(0,3)) +
  RotatedAxis()

publish_source_plot('psoriasis.pbmc.m2.dotplot')

sobj |>
  DotPlot2d('TRPM2', Status, manual_coarse) +
  labs(x = 'Group', y = 'Cell type') +
  theme_jpub() +
  scale_radius(range = c(0,3)) +
  RotatedAxis()

publish_source_plot('psoriasis.pbmc.m2.dotplot2')

sobj_mo |>
  DotPlot(c('TRPM2','CD14','FCGR3A'), cols = 'RdBu')

sobj_mo |> DimPlot(cols = DiscretePalette(36))

m2l_bc <- sobj_mo |>
  filter(seurat_clusters %in% c(5,11,15)) |>
  colnames()

sobj <- sobj |>
  mutate(manual_fine = case_when(.cell %in% m2l_bc ~ 'CD16+ Monocyte',
                                 manual_main == 'Monocytes' ~ 'CD14+ Monocyte',
                                 .default = manual_coarse))

sobj |>
  filter(manual_main != 'pDC') |>
  DotPlot2d('TRPM2', Status, manual_fine) +
  labs(x = 'Group', y = 'Cell type') +
  theme_jpub() +
  scale_radius(range = c(0,3)) +
  RotatedAxis()

publish_source_plot('psoriasis.pbmc.cd14mo.m2.dotplot', width = 60)

pso_pbmc_m2 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/psoriasis.pbmc.cd14mo.m2.dotplot.csv')

pso_pbmc_m2 |>
  mutate(group.y = fct_reorder(group.y, pct.exp, max)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 in PBMC') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('psoriasis.pbmc.cd14mo.m2.dotplot', width = 60)

## full rds ---------
sobj |>
  write_rds('mission/SLE_TRPM2_MfMo/psoria_pbmc.rds')

sobj_mo <- sobj |>
  filter(manual_main == 'Monocytes')

## mono ---------------
sobj_mo |>
  DotPlot('TRPM2', cols = 'RdBu') |>
  pluck('data') |>
  BubblePlot() +
  theme_jpub()

publish_source_plot('psoria.pbmc.mono.m2.dotplot')

sobj_mo |> DotPlot2d('TRPM2', Disease, seurat_clusters)

sobj_mo |> FeaturePlot('TRPM2', split.by = 'Disease', order = T,
                       cols = c('lightgrey','red')) &
  theme_jpub(theme_classic) & NoLegend()

publish_pdf('psoria.pbmc.mono.m2.featureplot.pdf', width = 100)

sobj_mo <- sobj_mo |>
  quick_process_seurat(skip_norm = T)

sobj_mo |> DotPlot(modc_marker, cols = 'RdBu', cluster.idents = T) +
  RotatedAxis()

sobj_mo |>
  FindMarkers(group.by = 'Status', ident.1 = 'Only psoriasis',
              ident.2 = 'Healthy', features = 'TRPM2')

sobj_mo |>
  FindMarkersAcrossVar(group.by = 'Disease', ident.1 = 'Psoriasis',
                       split.by = 'seurat_clusters', features = 'TRPM2')

### rds --------------
sobj_mo |>
  write_rds('mission/SLE_TRPM2_MfMo/psoria_pbmc_mono.rds')

sobj_mo <-
  read_rds('mission/SLE_TRPM2_MfMo/psoria_pbmc_mono.rds')

### mono umap ----------
sobj_mo |>
  DimPlot(cols = DiscretePalette(36), label = T, label.box = T, repel = T,
          label.size = 2) +
  theme_jpub(theme_classic)

publish_pdf('psoria.pbmc.mono.umap.pdf', width = 70)

c169_mrk <- sobj_mo |>
  FindMarkers(ident.1 = c(1,6,9), ident.2 = c(2,4,7,3,10,8))

c169_mrk |> filter(avg_log2FC > 0) |> head()

sobj_mo |> DotPlot(cc.genes.updated.2019, cols = 'RdBu', cluster.idents = T)

sobj_mo <- sobj_mo |>
  mutate(manual_fine = case_when(seurat_clusters == 13 ~ 'cDC',
                                 seurat_clusters %in% c(15,11,5) ~ 'ncMono',
                                 seurat_clusters %in% c(1,6,9,12) ~ 'cycling cMono',
                                 .default = 'cMono'))

sobj_mo |> DimPlot(group.by = 'manual_fine', cols = 'Paired') +
  ggtitle('PBMC monocytes')

sobj_mo |>
  bill.violin('TRPM2', group.by = Status) +
  labs(y = 'Normalized expression', title = 'Total monocytes')

sobj_mo |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'Status',
                       ident.1 = 'Psoriasis', features = 'TRPM2')

sobj_mo |>
  filter(manual_fine != 'cDC') |>
  DotPlot2d('TRPM2', Disease, manual_fine) +
  labs(x = 'Group', y = 'Monocyte subset')

sobj_mo |>
  DotPlot('TRPM2', cols = 'RdBu', cluster.idents = T)

sobj_mo <- sobj_mo |>
  AddModuleScore(features = list(inflam15$gene), name = 'inflam')

## inflam score -------------
inflam_m2 <- sobj_mo |>
  DotPlot(c('inflam1', 'TRPM2', 'CCR2', 'HMGB1')) |>
  pluck('data') |>
  as_tibble()

id2type <- sobj_mo |>
  distinct(seurat_clusters, manual_fine)

inflam_m2 |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  left_join(id2type, join_by(id == seurat_clusters)) |>
  filter(id != 13) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = manual_fine)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2, output.type = 'tex') +
  labs(title = 'Monocyte clusters in psoriasis PBMC', color = 'Subtype',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_jpub()

publish_source_plot('psoria.PBMC.m2.inflam.corr', width = 70)

sobj_mo <- sobj_mo |>
  filter(manual_fine != 'cDC') 

sobj_mo |>
  bill.violin(inflam15$gene, manual_fine, facet.ncol = 5) +
  labs(x = 'Cell type', fill = 'Cell type', y = 'Normalized expression') +
  theme_jpub(theme_classic) +
  theme(axis.text.x = element_blank())

publish_pdf('psoria.pbmc.mo.inflam.gene.violin.pdf', width = 90)

sobj_mo |>
  mutate(manual_fine = ifelse(manual_fine=='ncMono', 'TRPM2-lo Mo', 'TRPM2-hi Mo')) |>
  bill.violin(inflam15$gene, manual_fine, facet.ncol = 5) +
  labs(x = 'Cell type', fill = 'Cell type', y = 'Normalized expression') +
  theme_jpub(theme_classic) +
  theme(axis.text.x = element_blank())

publish_pdf('psoria.pbmc.m2hlmo.inflam.gene.violin.pdf', width = 90)

## pathway enrichment --------------
m2hvl_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'Disease', group.by = 'manual_fine',
                       ident.1 = 'ncMono') |>
  mutate(avg_log2FC = -avg_log2FC)

pso_m2hvl_gsego <- m2hvl_deg |>
  filter(cluster == 'Psoriasis', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(T) |>
  clusterProfiler::gseGO(OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL',
                         eps = 0)

pso_m2hvl_gsego <- pso_m2hvl_gsego |>
  clusterProfiler::simplify()

pso_m2hvl_gsego@result |>
  filter(NES > 0) |>
  plot_enrichment(metric = NES) +
  theme_jpub(theme_classic) +
  labs(title = 'GO BP pathway enrichment of psoriasis monocyte\nTRPM2-hi vs TRPM2-lo')

publish_source_plot('psoria.pbmc.m2hvl.gogse', width = 60)

## PBMC umap ------
pbmc_tibble <-
  read_csv('mission/SLE_TRPM2_MfMo/results/psoria_pbmc_meta.csv.csv')

pbmc_tibble |>
  ggplot(aes(x = umap_1, y = umap_2, fill = manual_main)) +
  geom_bin2d(bins = 512) +
  scale_fill_brewer(palette = 'Paired') +
  labs(fill = 'Cell type', title = 'Psoriasis patient PBMC') +
  theme_jpub(theme_classic)

publish_pdf('psoria.pbmc.umap.pdf', width = 70)

## cell fraction --------
mosub_bysample <- sobj_mo |>
  calc_frac_conf_on_grouped_count(Subject, manual_fine)

mosub_bysample |>
  mutate(group = ifelse(str_detect(Subject, 'HC'), 'HC', 'Psoriasis')) |>
  ggplot(aes(group, fraction*100)) +
  geom_boxplot() +
  facet_grid(~manual_fine) +
  stat_compare_means(method = 't.test', label = 'p.signif') +
  labs(y = '% in total monocytes')

mosub_bysample |>
  mutate(group = str_extract(Subject, 'HC|PSA|PSO')) |>
  ggplot(aes(group, fraction*100)) +
  geom_boxplot() +
  facet_grid(~manual_fine) +
  stat_compare_means(method = 't.test', size = 2,
                     comparisons = list(c('HC','PSA'),c('HC','PSO'))) +
  labs(y = '% in total monocytes') +
  theme_jpub()

publish_source_plot('psoria.pbmc.monosub.frac', width = 70)

mosub_bysample <- sobj_mo |>
  calc_frac_conf_on_grouped_count(Subject, seurat_clusters)

mosub_bysample |>
  mutate(group = ifelse(str_detect(Subject, 'HC'), 'HC', 'Psoriasis')) |>
  filter(seurat_clusters == 5) |>
  ggplot(aes(group, fraction)) +
  geom_boxplot() +
  stat_compare_means(method = 't.test')

pbmc_bysample <- pbmc_tibble |>
  calc_frac_conf_on_grouped_count(Subject, CellType)

pbmc_bysample |>
  mutate(group = ifelse(str_detect(Subject, 'HC'), 'HC', 'Psoriasis'),
         CellType = case_when(str_detect(CellType, 'CD14') ~ 'TRPM2-hi Mono',
                              str_detect(CellType, 'CD16') ~ 'TRPM2-lo Mono',
                              .default = NA)) |>
  na.omit()|>
  ggplot(aes(group, fraction*100)) +
  geom_boxplot() +
  stat_compare_means(method = 't.test') +
  facet_grid(~CellType) +
  labs(y = '% in total PBMC')

pbmc_bysample |>
  mutate(group = str_extract(Subject, 'HC|PSA|PSO')) |>
  filter(str_detect(CellType, 'ono')) |>
  ggplot(aes(group, fraction)) +
  geom_boxplot() +
  stat_compare_means(method = 't.test', ref.group = 'HC', label = 'p.signif') +
  facet_grid(~CellType)
